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Abstract 

Modern  theories  of  gravity  prediet  ripples  in  space  and  time,  which  are  known  as 
gravitational  waves.  Very  large  interacting  masses,  sueh  as  binary  systems  of  neutron  stars  or 
black  holes,  are  predieted  to  generate  gravitational  waves  that  may  be  intense  enough  to  be 
deteeted.  This  study  presents  a  response  analysis  for  the  proposed  gravitational  wave  experiment 
known  as  the  Laser  Interferometer  Space  Antenna  (LISA).  The  analysis  focuses  on  the 
gravitational  wave  distortion  patterns,  or  polarizations,  predicted  by  the  Brans-Dicke  Scalar- 
Tensor  Theory  of  Gravity  (BD).  Several  plausible  metric  theories  of  gravity,  sueh  as  BD,  predict 
gravitational  wave  polarizations  that  differ  from  the  traditionally  expeeted  polarizations 
predicted  by  General  Relativity  (GR).  LISA,  a  cooperative  venture  between  the  National 
Aeronauties  and  Spaee  Administration  (NASA)  and  the  European  Space  Ageney  (ESA),  is  a 
proposed  space-based  gravitational  wave  experiment  with  the  potential  to  distinguish  between 
polarizations  of  gravitational  waves.  EISA  is  scheduled  for  launch  in  2012;  so  much  of  the 
eurrent  work  for  EISA  involves  the  development  and  design  of  systems  and  procedures  that  will 
allow  EISA  to  attain  its  desired  performance. 

Eor  a  gravitational  wave  experiment,  the  projected  sensitivity  is  determined  by  an 
instrument  response  function,  which  relates  the  disturbanee  generated  by  the  gravitational  wave 
to  the  disturbance  that  can  be  measured  by  the  instrument.  This  study  compares  the  response 
funetions  for  BD  polarizations  to  the  published  response  functions  for  GR.  Response  functions 
corresponding  to  several  different  methods  of  data  eollection  were  analyzed.  The  analysis 
revealed  that  eaeh  response  function  depends  strongly  on  the  polarization,  frequency,  and 
propagation  direction  of  the  gravitational  wave. 

To  numerieally  evaluate  and  eompare  the  response  funetions  two  programs  were 
generated  in  the  MATLAB®  programming  language.  The  first  program  computes  the  response 
functions  averaged  over  all  souree  direetions  for  eaeh  gravitational  wave  polarization  and  data 
colleetion  method.  The  seeond  program  caleulates  the  instrument  response  at  eaeh  point  in 
EISA’s  orbit  to  gravitational  waves  originating  from  the  galaetic  eenter.  The  results  from  the 
programs  generated  in  Matlab  identify  patterns  in  the  response  funetions  that  distinguish  between 
gravitational  waves  of  BD  and  GR  polarizations.  These  results  identify  features  of  the  response 
functions  that  may  assist  future  researchers  in  determining  the  optimal  data  colleetion  method  to 
deteet  gravitational  waves  of  a  given  polarization,  frequency,  and  source  direction. 
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Introduction 


The  discovery  and  analysis  of  gravitational  waves  promises  to  usher  in  a  new  age  of 
astronomy.  Nearly  all  astronomical  data  are  gathered  from  electromagnetic  radiation  emitted  by 
accelerated  electrical  charges  in  distant  stars.  The  detection  and  subsequent  analysis  of  radiation 
emitted  by  accelerated  mass  would  provide  a  new  probe  of  the  universe  and  its  structure. 

The  apparent  parallels  between  the  gravitational  interaction  and  the  electromagnetic  force 
have  long  been  recognized.  The  most  obvious  parallel  is  the  similarity  between  Isaac  Newton’s 
Law  of  Universal  Gravitation  and  Charles  Coulomb’s  Electrostatic  Force  Law: 


^  Gin.irij 

F  = - ^ 

s  2 


F  = 


kq,q^ 


(1) 

(2) 


Both  equations  involve  a  constant  (G  or  k)  multiplied  by  a  fundamental  quantity  (mass 
and  charge)  of  two  objects,  and  divided  by  the  distance  separating  the  objects  squared.  Until  the 
20*  century,  no  one  considered  searching  for  a  gravitational  radiation  similar  to  electromagnetic 
radiation  but  generated  by  the  acceleration  of  mass.  Because  Newton’s  theory  of  gravity  fails  to 
support  the  wave  equation  necessary  for  the  existence  of  gravitational  radiation,  a  search  for 
gravitational  radiation  was  not  merited  until  the  general  theory  of  relativity  was  published  in 
1915  [1]. 

In  Newton’s  theory,  gravity  is  considered  to  be  a  force  that  is  experienced 
instantaneously,  that  is  to  say  that  if  the  Sun  were  to  suddenly  vanish,  the  earth  and  all  the  other 
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planets  would  immediately  depart  their  orbits  and  head  for  deep  space  with  the  velocity  they  had 
the  moment  the  sun  vanished.  Einstein  departed  from  the  instantaneous  gravitational  effect 
intrinsic  to  Newton’s  theory  by  assuming  that  no  information,  including  the  gravitational 
interaction,  could  propagate  faster  than  the  speed  of  light  in  vacuum,  which  he  assumed  as 
invariant  [2], 

Einstein  thoroughly  explored  the  consequences  of  assuming  the  speed  of  light  to  be  a 
constant,  and  published  the  results  of  his  investigation  in  1905  in  his  Special  Theory  of  Relativity 
[1].  Special  relativity  allows  one  to  find  the  relationship  between  measurements  made  in  two 
inertial  reference  frames.  To  illustrate  the  concept  of  special  relativity  (See  Eig.  1),  consider  a 
situation  where  Bob  stands  on  a  flatbed  railcar  with  a  timepiece  that  consists  of  light  that 
bounces  between  two  mirrors. 


Figure  1 :  In  reference  frame  A,  Bob  observes  a  timepiece  consisting  of  a  photon  bouncing 
between  two  mirrors. 


Eight  has  the  special  property  that  it  travels  through  space  at  a  constant  rate  (c)  that  does 


not  depend  on  the  motion  of  the  source  or  receiver.  As  the  light  completes  one  round  trip.  Bob 
(in  reference  frame  A)  will  observe  the  light  traveling  a  distance  2d  in  the  time  t. 
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where  ^  .  Meanwhile,  in  reference  frame  B,  Alice  stands  on  the  ground  and  watches  the 

flatbed  car  travel  by  her  at  a  velocity  v  (See  Fig.  2). 


Figure  2:  In  reference  frame  B,  Alice  observes  the  path  taken  by  a  photon  traveling  between  two 
mirrors  that  are  moving  relative  to  her  position. 


Alice  observes  the  light  traveling  a  distance  greater  than  2d  as  it  makes  a  round  trip. 
Consequently,  she  notices  more  time  has  elapsed  between  the  same  two  events  that  Bob 
witnessed  in  Reference  Frame  A.  Though  the  amount  of  time  observed  between  two  events  is 
dependent  upon  the  reference  frame,  Einstein  concluded  that  the  right  mixture  of  space  and  time, 
referred  to  as  the  spacetime  interval,  does  not  depend  upon  the  reference  frame.  The  spacetime 
interval  and  the  principle  that  there  is  no  preferred  frame  or  observer  are  the  essence  of  special 
relativity.  The  invariant  interval  is  given  by  the  equation: 

Ay^  =  Ac^  +  Ap^  +  Az^  -  A(cO^  ■  (3) 

Therefore,  light  in  all  reference  frames  propagates  between  spacetime  points  separated  by 

2 

a  null  interval  {As  =  0).  Special  relativity  is  only  valid  for  reference  frames  that  are  moving  at 
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constant  velocities  and  does  not  aeeount  for  gravitational  effeets.  In  order  to  apply  the  preeepts 
of  relativity  to  a  gravitational  context,  Einstein  used  what  is  called  the  equivalence  principle  [2]. 
To  begin,  Einstein  surmised  that  the  effeets  of  gravity  were  no  different  than  the  effects 
experienced  in  an  accelerated  reference  frame.  Eor  example,  consider  the  two  reference  frames 
illustrated  in  Fig.  3. 


Figure  3:  A  rocket  ship  blasts  upwards  so  an  observer  would  witness  laser  light  inside  the 
spacecraft  following  a  curved  line  (A).  The  acceleration  due  to  the  gravitation  of  the  earth  causes 
light  to  bend  in  the  same  manner  as  seen  on  an  accelerating  spacecraft  (B). 

The  two  reference  frames  consist  of  identical  rooms;  one  room  located  on  a  spacecraft 
accelerating  upwards  at  9.8  m/s  and  the  other  room  located  at  the  surface  of  the  earth  where  the 
gravity  causes  a  free  fall  acceleration  of  9.8  m/s  .  In  both  rooms  a  laser  at  the  left  end  of  the 
room  shoots  a  beam  that  hits  the  right  side  of  the  room.  Because  the  rocket  in  reference  frame  A 
accelerates  up  during  the  time  the  laser  light  travels  across  the  room,  the  path  of  the  light  appears 
to  bend  so  that  it  impacts  the  wall  on  the  right  side  of  the  room  at  a  point  lower  then  the  initial 
height  of  the  laser.  The  equivalence  principle  asserts  that  a  gravitational  field  will  have  the  same 
effect  on  the  laser  beam.  Einstein  realized  that  the  path  taken  by  the  laser  could  be  described  by 
a  geodesic  curve  in  spacetime.  He  was  then  able  to  create  the  general  theory  of  relativity  by 
developing  equations  that  describe  how  matter  influences  the  geometry  of  spacetime  while  the 
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geometry  of  spacetime  influences  the  motion  of  matter  [1],  The  equations  Einstein  developed, 
known  as  field  equations,  are  analogous  to  Maxwell’s  equations  for  electromagnetism,  and  they 
predict  the  wave  equation  that  Newton’s  theory  did  not  support. 

Because  Einstein’s  theory  describes  gravity  in  terms  of  the  geometry  of  spacetime,  it  is 
referred  to  as  a  metric  theory  of  gravity.  After  studying  Einstein’s  approach  to  gravity,  several 
other  scientists  have  developed  alternative  metric  theories  that  describe  the  connection  between 
matter  and  spacetime  curvature  differently.  In  the  1960’s  Carl  Brans  and  Robert  Dicke 
published  a  scalar-tensor  theory  of  gravity  [3].  The  Brans-Dicke  theory  (BD)  was  developed  to 
incorporate  Mach’s  principle  into  Einstein’s  theory  of  gravity.  Mach’s  principle  asserts  that  the 
inertial  properties  of  matter  are  dependent  upon  the  total  distribution  of  matter  in  the  universe 
[3].  The  essential  result  of  BD  is  that  the  gravitational  “constant”  G  is  not  constant,  but 
dependent  upon  the  totality  of  matter  in  the  universe  and  can  vary  from  place  to  place  and  time  to 
time  [3].  In  general,  BD  and  other  alternative  theories  agree  with  the  results  of  GR  when  the 
gravitational  field  is  relatively  weak,  such  as  within  the  solar  system.  However,  as  the 
gravitational  field  becomes  stronger,  the  mathematical  description  of  spacetime  curvature 
becomes  increasingly  complicated,  and  the  predictions  of  alternative  theories  diverge  from  the 
predictions  of  GR  [4].  Gravitational  waves  with  sufficient  amplitude  for  detection  are  only 
created  within  the  strongest  gravitational  fields.  As  a  gravitational  wave  leaves  the  strong  field, 
the  persisting  polarization  and  frequency  of  the  wave  contains  information  on  the  nature  of  the 
strong  field  produced  by  the  source.  Because  it  is  not  possible  to  send  a  gravitational  experiment 
to  the  surface  of  a  distant  neutron  star  or  black  hole  where  the  gravitational  field  is  recognized  as 
strong,  the  properties  of  the  gravitational  interaction  in  the  strong  field  are  not  well  known.  The 
detection  and  analysis  of  gravitational  waves  may  be  the  first  test  to  validate  theories  that  differ 
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from  GR  only  in  the  limit  of  the  strong  field  [5].  Gravitational  wave  deteetion  stands  as  the  long 
awaited  window  into  the  nature  of  the  strongest  gravitational  fields. 

Theoretieally,  gravitational  waves  originate  with  any  aeeeleration  of  mass;  however,  due 
to  the  weak  nature  of  the  gravitational  interaetion,  only  the  most  intense  aeeelerations  of  the 
densest  matter  will  result  in  gravitational  waves  that  may  be  measurable  [6].  Binary  systems  of 
neutron  stars  or  blaek  holes  are  likely  to  produee  gravitational  waves  of  suffieient  strength  to 
measure  [7].  Gravitational  waves  aet  to  streteh  or  shrink  the  spaee  that  they  travel  through  in 
predietable  patterns.  Consequently,  objeets  in  the  path  of  a  gravitational  wave  will  experienee  a 
eompressive  or  tensile  strain,  and  the  distanees  between  objeets  will  either  deerease  or  inerease. 
Strain  is  measured  by  the  ehange  in  length  over  the  original  length  (AL/L).  The  distortion 
pattern,  or  polarization,  of  the  gravitational  wave  determines  the  direetion  an  objeet  will  be 
strained,  and  similarly  the  direetion  objeets  will  move  in  relation  to  one  another.  GR  prediets 
two  polarizations,  whieh  are  known  as  the  plus  and  eross  polarizations  [8].  In  addition  to  the  two 
polarizations  predieted  by  GR,  BD  prediets  the  existenee  of  a  third  polarization,  designated  here 
as  the  sealar  polarization  [6].  All  three  polarizations  eause  spatial  strains  that  are  transverse 
relative  to  the  direetion  of  propagation.  Figure  4  depiets  how  a  gravitational  wave  with  a  plus 
polarization  would  affeet  a  group  of  masses  arranged  in  a  eirele.  The  distortion  pattern  is 
mapped  out  over  one  period  in  the  gravitational  wave. 
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Plus  Polarization 

Propagation  Direction  Out  of  the  Page 


0  T/4  T/2  3T/4  T 


Figure  4:  The  strain  pattern  exhibited  by  the  plus  polarization  in  one  cycle  of  a  gravitational  wave 
illustrated  in  a  series  of  snapshots.  The  numbers  below  each  snapshot  indicate  the  progression  of 
the  gravitational  wave  through  a  single  cycle  -  zero  indicates  the  beginning  of  a  cycle  and  T 
indicates  the  end  of  the  cycle.  [6] 


The  cross  polarization  behaves  similarly  although  the  strain  pattern  is  rotated  by  45  degrees 
about  the  direction  of  propagation. 


Cross  Polarization 

Propagation  Direction  Out  of  the  Page 


0  T/4  T/2  3T/4  T 


Figure  5:  The  strain  pattern  exhibited  by  the  cross  polarization  through  one  cycle.  [6] 


Finally,  the  scalar  polarization  acts  as  a  breathing  mode  causing  a  tensile  strain  in  all  transverse 
directions  followed  in  time  by  a  compressive  strain. 
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Scalar  Polarization 

Propagation  Direction  Out  of  the  Page  0 


0  T/4  T/2  3T/4  T 


Figure  6:  The  strain  pattern  exhibited  by  the  scalar  polarization  through  one  cycle.  [6] 

The  illustrations  of  the  three  polarizations  are  helpful  to  visualize  the  strain  pattern  of  a 
gravitational  wave;  however,  they  in  no  way  indieate  the  amplitude  of  the  wave.  Known  systems 
of  binary  neutron  stars  are  predicted  to  produce  gravitational  waves  that  will  reach  earth  with  a 
strain  amplitude,  a  fractional  change  in  length  over  length,  on  the  order  of  10'  [7],  which  is  the 
approximate  design  sensitivity  for  LISA.  Such  a  gravitational  wave  will  change  the  length  of  a 
thousand  kilometer  ruler  by  about  the  diameter  of  a  proton.  Two  approaches  have  been  explored 
to  detect  such  a  small  strain.  The  first  involves  a  solid  aluminum  bar  designed  to  “ring”  when  a 
gravitational  wave  with  a  particular  frequency  travels  through  it.  If  the  aluminum  bar’s 
fundamental  mode  of  vibration  has  a  frequency  that  matches  that  of  the  gravitational  wave,  then 
each  passing  phase  of  the  gravitational  wave  will  add  to  the  vibrational  energy  of  the  bar  in  a 
resonance  effect.  The  resulting  vibrational  strain  in  the  bar  could,  in  theory,  be  several  orders  of 
magnitude  greater  then  the  amplitude  of  the  incident  wave.  However,  no  bar  detector  has 
verifiably  detected  a  gravitational  wave  [9]. 
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The  second  approach  uses  lasers  to  measure  changes  in  the  distance  between  points  in 
space  that  are  marked  off  by  proof  masses.  The  proof  masses  can  be  thought  of  as  points  on  the 
circles  used  to  illustrate  the  different  polarizations.  For  ground  based  experiments  the  proof 
masses  define  points  at  the  end  of  arms  that  are  mutually  perpendicular  to  one  another 
(See  Fig.  7). 


Ground  Based  Interferometer 


Figure  7:  A  depiction  of  an  earth  based  gravitational  wave  interferometer.  As  the  gravitational 
wave  warps  the  space  containing  the  interferometer  the  lengths  of  the  respective  arms  will  shrink 
or  grow  by  an  amount  dictated  by  the  amplitude,  polarization,  and  propagation  direction  of  the 
gravitational  wave.  [9] 


A  laser  beam  is  sent  through  a  beam  splitter  and  directed  down  two  paths  that  are  equal  in 
length  (Lx  and  Ly).  This  setup  will  be  referred  to  in  this  report  as  the  two-arm  Michelson  mode 
of  detection.  The  beams  reflect  off  the  proof  masses  and  recombine  at  the  beam  splitter  where 
they  are  directed  to  a  detector.  The  intensity  of  the  light  measured  at  the  detector  is  related  to  the 
difference  in  the  lengths  of  each  arm,  |Lx  -  Ly|.  If  both  paths  are  exactly  equal  in  length  then  the 
beams  will  destructively  interfere  and  the  detector  will  measure  an  intensity  minimum.  If  one 
path  is  longer  then  the  other  by  half  of  the  wavelength  of  the  laser  light,  then  the  beams  will 
constructively  interfere  and  the  detector  will  measure  an  intensity  maximum.  The  shift  from  an 
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intensity  minimum  to  maximum  is  known  as  one-half  of  an  interferenee  fringe  shift.  The  ground 
based  interferometer  experiments  are  limited  by  the  length  of  the  arms  (approximately  4  km)  and 
by  seismie  noise  which  disturbs  the  instrument  causing  false  signals  [9].  A  similar 
interferometer  built  in  space  could  have  arm  lengths  nearly  500  times  the  diameter  of  the  earth 
and  would  not  be  subject  to  seismic  noise  [7]. 

The  design  and  construction  of  a  space-home  interferometer  known  as  the  Laser 
Interferometer  Space  Antenna  (LISA)  is  currently  being  led  by  the  National  Aeronautics  and 
Space  Administration’s  Goddard  Space  Flight  Center.  In  the  LISA  design,  spacecraft  will 
occupy  the  vertices  of  an  equilateral  triangle  5  million  kilometers  on  a  side.  The  significant 
increase  in  path  length  reduces  the  fractional  change  necessary  to  give  a  fringe  shift.  As  a  result, 
the  effective  detection  band  is  shifted  to  lower  frequencies.  LISA  is  designed  for  the  10'"^  to  10*’ 
Hertz  band.  The  triangle  will  follow  the  earth  in  its  orbit  and  will  rotate  clockwise  as  it  revolves 
around  the  sun  [7]  (See  Fig.  8). 


A  B 

Figure  8:  These  pietures  illustrate  the  LISA  configuration.  Picture  A  illustrates  how  each  spacecraft 
follows  its  own  orbit  about  the  sun  so  that  the  three  spacecraft  will  maintain  relative  positions  at  the 
vertices  of  an  equilateral  triangle  throughout  the  course  of  a  complete  revolution  about  the  sun. 
Picture  B  shows  the  location  of  the  LISA  configuration  relative  to  the  earth’s  orbit  [10]. 
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Each  spacecraft  contains  a  laser  and  proof  masses  so  the  eonstellation  of  three  spaeeeraft 
ean  operate  as  one  instrument  on  the  same  principle  as  the  ground  based  experiments.  Beeause 
eaeh  spaeeeraft  in  the  LISA  setup  contains  an  interferometer,  the  data  gathered  from  the  fringe 
shifts  measured  at  the  detectors  can  be  related  to  several  different  combinations  of  the  arm 
lengths  [11]. 


Figure  9:  The  LISA  reference  frame  defined  in  this  study.  LI,  L2,  and  L3  represent  the  lengths  of 
each  arm. 

For  example,  in  referenee  to  Fig.  9,  spaeeeraft  A  will  collect  data  corresponding  to  |L1- 
L3|,  while  spaeeeraft  B  will  colleet  data  corresponding  to  |L1-L2|,  and  so  on.  Data  transmitted  to 
earth  can  be  reproeessed  to  yield  other  eombinations  including  permutations  of  |L1+L2-2L3|  and 
|L1+L2+L3|  [11].  The  different  kinds  of  eombinations  will  be  referred  to  in  this  paper  as  modes 
of  data  colleetion.  Four  modes  of  data  collection  are  examined  in  this  study.  The  first  mode 
evaluates  the  response  of  a  single  arm,  the  seeond  computes  the  LISA  equivalent  to  the  two-arm 
Michelson  mode  that  is  used  in  earth  based  detectors.  The  two  arm  Michelson  is  defined  by  |L1- 
L2|.  The  third  mode  is  a  result  of  eombining  the  results  of  two  independent  two  arm  Michelson 
responses.  It  is  referred  to  as  the  three-arm  Michelson  mode  and  is  defined  by  |L1+L2-2L3|. 

The  final  mode  utilizes  the  LISA  setup  as  a  Sagnae  interferometer  so  it  is  appropriately  titled  the 
Sagnae  mode  and  is  designated  as  |L1+L2+L3|. 
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Determining  the  sensitivity  of  eaeh  mode  of  data  eolleetion  to  gravitational  waves  of 
various  frequeneies,  source  directions  and  polarizations  is  an  essential  part  of  the  LISA  program. 
Several  groups  have  published  studies  of  LISA’ s  sensitivity  to  the  gravitational  wave 
polarizations  predicted  by  GR  [7,  11,  12],  Other  groups  have  reported  that  no  experiment  to  date 
provides  data  to  exclude  scalar-tensor  theories  like  BD.  The  following  report  is  the  first  to 
explore  the  response  of  LISA  to  the  additional  gravitational  wave  polarization  predicted  by  BD. 
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Methods 


The  sensitivity  of  LISA  ean  be  ealculated  through  a  eombination  of  two  funetions;  a 
noise  funetion  and  a  response  funetion.  The  noise  function  is  a  measure  of  the  equivalent  strain 
induced  by  instrument  noise.  Self  gravity  noise,  residual  gas  damping  and  minor  fluctuations  in 
the  laser  amplitude  and  frequency  can  simulate  strain  that  could  be  mistaken  for  a  passing 
gravitational  wave  [11].  The  response  function  relates  the  strain  produced  by  a  gravitational 
wave  to  the  time  delay  experienced  by  the  beams  propagating  between  pairs  of  spacecraft  in  the 
LISA  system.  The  two  functions  are  combined  to  form  a  plot  used  to  determine  the  signal  to 
noise  expected  for  LISA. 

The  noise  function  is  dependent  only  upon  instrument  parameters  and  the  detector 
environment.  Neither  the  polarization  nor  the  frequency  nor  any  other  aspect  of  a  gravitational 
wave  shows  up  in  the  noise  function.  Consequently,  every  study  of  LISA’ s  sensitivity, 
regardless  of  the  polarization  or  source  direction  examined  should  contain  the  same  noise 
function.  Therefore,  the  main  focus  of  this  study  is  to  develop  and  analyze  the  response  function 
of  LISA  for  the  scalar  polarization  predicted  in  BD. 

Systems  engineers  are  often  concerned  with  developing  a  mathematical  formula  that 
relates  the  input  to  a  system  to  the  expected  output.  The  formula  sought  is  often  referred  to  as  a 
transfer  function.  A  response  function  is  very  similar  to  a  transfer  function  in  the  way  it  relates 
the  gravitational  wave,  which  may  be  considered  the  input,  to  the  spatial  strain  experienced  by 
the  laser  light  traveling  between  two  spacecraft.  To  better  understand  the  response  function, 
consider  an  analogous  situation  involving  corks  floating  on  a  pond. 
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Figure  10:  The  distance  between  two  corks  changes  as  they  bob  up  and  down  on  waves  across  a 
pond.  The  green  line  represents  a  hypothetical  rubber  band  used  to  measure  the  changes  in 
distance  between  the  corks. 

A  lightweight  rubber  band  stretched  between  the  corks  stretches  as  a  wave  on  the  pond 
causes  the  corks  to  move  up  and  down.  If  details  of  the  wave  moving  the  corks  are  known,  i.e. 
the  amplitude,  frequency,  direction  and  speed,  then  it  would  be  possible  to  generate  a 
mathematical  model  to  calculate  how  much  the  rubber  band  will  stretch  when  disturbed  by  a 
passing  wave.  The  response  function  is  the  backbone  of  such  a  model  because  it  is  able  to 
incorporate  the  effects  of  different  frequencies  and  source  directions  of  the  incident  wave.  A 
wave  originating  from  the  right  or  left  of  the  corks  depicted  in  Fig.  10  would  produce  a  response 
that  varies  in  time.  In  the  case  where  the  waves  were  originating  from  above  or  below  the  page 
(Fig.l  1),  both  corks  could  simultaneously  be  at  the  crest  or  the  trough  of  the  wave  and  the  rubber 


band  would  not  stretch  at  ah. 
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Figure  1 1 :  The  direction  of  the  incident  waves  affects  how  the  distance  between  the  corks  changes. 

In  this  picture,  the  waves  are  rolling  into  the  page  while  both  corks  are  bobbing  up  and  down  in 
unison. 

The  situation  in  which  the  rubber  band  does  not  stretch  is  termed  a  null  response.  The 
frequency  of  the  wave  may  also  produce  a  null  response;  in  fact,  as  the  frequency  of  the  wave 
increases,  the  wavelength  will  become  shorter  until  it  matches  the  distance  between  the  corks. 
At  that  point  both  corks  would  reside  on  different  crests  and  the  rubber  band  would  not  be 
stretched  (Fig.  12).  The  response  function  should  indicate  a  series  of  null  responses  as  the 
frequency  continues  to  increase  and  the  distance  between  the  corks  becomes  equal  to  integer 
multiples  of  the  wavelength. 


Figure  12:  The  frequency  of  incident  waves  also  affects  the  response  function.  The  distance 
between  corks  does  not  change  because  the  corks  float  on  the  same  location  of  different  wave 
crests. 
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Before  an  equation  ean  be  developed  to  deseribe  the  response  funetion  of  LISA,  the 
gravitational  wave  must  be  deseribed  mathematieally.  A  gravitational  wave  ean  be 
approximated  as  a  perturbation  of  spaeetime,  so  the  deseription  of  a  gravitational  wave  must 
begin  with  the  deseription  of  spaeetime.  A  relatively  small  region  of  spaee  and  time  ean  be 
eonsidered  flat  to  a  first  order  approximation.  The  spaeetime  interval  for  a  flat  region  may  be 
ealeulated  using  the  metrie  of  speeial  relativity  (See  Eq.  3)  written  here  in  differential  form; 

ds^  =  dx^  +  dy^  +  dz^  -  dt^ ,  c  =  l .  (4) 


The  spaeetime  interval  may  also  be  written  in  metrie  form; 
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The  flat  metrie 
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is  symbolized  as  where  p  and  v  are  indiees  that  range  in  value  from  0 


to  3  [8].  Therefore, 
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Using  Einstein’s  summation  notation,  the  summations  over  |u  and  v  are  implied  by  the  repeated 
indie  es  and  no  ^  is  used; 


ds^  =  Tj^^dx^dx'' .  (7) 

Equation  7  ean  then  be  modified  to  represent  a  gravitational  wave.  Eor  example,  the  spaeetime 
interval  for  a  perturbation  to  flat  spaeetime  by  a  transverse  gravitational  wave  propagating  along  the 
z  axis  is  given  by 


Here  the  perturbation  matrix  is: 
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Eor  the  sake  of  simplifying  the  notation,  eoeffieient  matriees  will  be  denoted  with  a  eapital  letter, 

i.e. 


Kv  =  ^  +  5) . 


(10) 
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The  amplitude,  A,  of  the  perturbation  must  necessarily  be  much  less  than  one,  which  is  consistent 
with  predicted  strain  amplitudes  of  gravitational  waves  (on  the  order  of  10'^').  In  matrix  form, 
the  coefficients  of  the  perturbation  matrices  for  the  cross,  plus,  and  scalar  polarizations  in  the 
source  frame  are  given  below. 
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Cross  Plus  Scalar 


GR  predicts  only  the  cross  and  plus  polarizations,  while  BD  predicts  all  three.  The  term 
‘source  frame’  is  in  reference  to  the  orientation  of  the  coordinate  axes  at  the  source,  where  the  z 
axis  is  reserved  for  the  direction  of  propagation.  It  is  important  to  note  that  LISA  must  be 
described  by  a  set  of  coordinate  axes  that  are  independent  of  those  in  the  source  frame 
(See  Fig.  13). 


Figure  13:  An  illustration  depicting  the  relationship  between  the  source,  heliocentric,  and  LISA 
reference  frames.  indicates  the  z  axis  of  the  source  frame  while  Zh,  and  Zl  indicate  the  z  axes 
of  the  heliocentric  and  LISA  frames  respectively. 
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Ignoring  the  orbital  motion  around  the  sun  and  assuming  LISA  remains  fixed  in  spaee, 
two  rotations  are  required  to  transform  from  the  source  frame  to  the  LISA  frame.  An  additional 
three  rotations  are  required  when  taking  the  orbital  motion  of  LISA  into  account  [7].  The 
process  of  calculating  the  transformations  is  outlined  in  appendix  A.  It  is  important  to  note  that 
the  rotations  only  affect  the  coefficient  matrix. 

With  the  description  of  the  gravitational  wave  in  the  LISA  coordinate  frame,  it  is  possible 
to  investigate  how  the  wave  influences  the  spacetime  traveled  by  the  laser  beams  passing 
between  the  spacecraft.  The  light  travels  along  a  null  geodesic,  which  is  a  path  in  spacetime 
where  the  spacetime  interval  is  equal  to  zero. 

0  =  +  H ^^Acos{cot  -  k  ■  f  +  S))dx^ dx'' ,  6  =  at  -  k  -  r  .  (11) 


Here  k  is  the  vector  associated  with  the  propagation  direction  and  wavelength  of  the 
gravitational  wave,  r  is  the  vector  associated  with  the  distance  traveled  by  the  light  along  an  arm 
of  LISA,  and  5  is  a  phase  constant.  Allowing  A  to  equal  one  will  produce  a  response  function 
that  indicates  the  fraction  of  the  wave’s  amplitude  that  is  measured  by  the  detector.  The  light 
travel  time  between  spacecraft,  t,  is  approximately  16  seconds.  For  waves  with  periods  shorter 
than  16  seconds,  the  space  is  alternately  stretched  and  compressed  as  the  light  pulse  travels 
between  spacecraft.  The  result  is  a  roll-off  in  the  response  for  frequencies  above  10"'  Hz. 
Separating  the  time  component  from  the  space  components  yields 


=  idnm  +  H„„Acos(cot  -k-r  +  S))dx’'dx"' . 


(12) 
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The  subscripts  n  and  m  are  indices  that  range  from  1  to  3,  and  include  only  the  spatial 
components  of  the  matrices.  Taking  the  square  root  of  both  sides  allows  for  an  integration  of  the 
light  travel  time.  Because  the  perturbations  to  spacetime  are  small,  the  binomial  expansion  may 
be  employed  to  rid  the  integral  of  the  radical. 


^0  ^0  1 

^dt  =  ^{t]^^+—H^^Acos{o}t-k-r  +  5))dx"dx'^  .  (13) 

The  integral  on  the  right  hand  side  is  evaluated  over  the  spacetime  path  of  a  light  pulse  traveling 
along  the  arm  of  LISA  being  examined.  Note  that  dx,  dy,  and  dz  can  be  defined  in  terms  of  the 
original  length  of  the  arm  when  the  LISA  frame  is  defined  (See  Fig.  14). 


Figure  14:  The  LISA  frame.  For  an  integration  along  hi  \  dx^  dt  cos(a),  dy  =  dt  sin(a),  and 
dz  =  0. 
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The  response  funetion  re  fleets  the  time  delay  for  light  traveling  an  arm  of  LISA  and  is  therefore 
given  by: 


t'—t 
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1 

jj  cos^(«)  +  //j2  eos(«)sin(«)  +  sin(«)cos(«)  +  sin^(«)]eos((yt  -k  ■r  +  S)dt .  (14) 

0 

The  remaining  integral  may  be  solved  analytically. 
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Results 


For  this  study,  MATLAB  was  used  to  produce  two  programs  to  numerically  evaluate  and 
compare  response  functions  generated  with  respect  to  the  three  different  polarizations  predicted 
by  BD.  The  first  program  computes  the  response  functions  for  each  polarization  and  data 
collection  mode  by  averaging  over  all  incident  source  directions.  The  second  program  calculates 
the  instrument  response  to  gravitational  waves  originating  from  the  galactic  center  as  LISA 
orbits  the  sun. 

To  validate  the  mathematics  and  programming  employed  to  calculate  the  response 
functions,  the  results  of  the  two-arm  Michelson  to  GR  polarizations  computed  by  the  averaging 
program  were  compared  to  those  published  by  Larson,  Hiscock,  and  Hellings  (LHH)  in  1999 
(See  Fig.  15).  Minor  discrepancies  between  the  plots  are  the  result  of  different  approaches  and 
mathematical  approximations.  The  LHH  study  involved  a  complicated  approach  to  calculate  the 
detector  response  in  terms  of  the  power  of  the  gravitational  wave,  while  the  approach  taken  in 
this  study  determines  the  response  in  terms  of  the  amplitude  of  the  gravitational  wave.  In 
addition,  the  LHH  approach  assumes  a  slightly  different  arm  length  reflected  as  a  small  shift 
along  the  log  (u)  axis.  The  logarithmic  plot  of  the  power  response  for  this  study  (illustrated  in 
Figure  15)  was  generated  by  taking  the  log  of  the  square  of  the  amplitude. 
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Averaged  P<wer  Response  Two  Arm  Michelson  Mode 
General  Relativity  Polarizations 


LHH  Response  Averaged 


Figure  15:  The  averaged  power  response  for  GR  polarizations  interaeting  with  LISA  in  the  Two 
Am  Miehelson  mode  found  in  this  study  (top)  is  eompared  to  the  power  response  of  LISA 
published  by  Shane  Larson,  William  Hiscoek  and  Ronald  Hellings  in  1999  (bottom)  [11].  The 
variable  u  represents  the  angular  ffequeney  {2nf). 
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Once  it  was  established  that  the  program  could  produce  reliable  results,  plots  of  the 
averaged  LISA  response  to  the  scalar  polarization  were  generated  for  comparison  to  the  averaged 
LISA  response  to  GR  polarizations.  The  response  functions  were  generated  for  each  mode  of 
data  collection,  which  includes  the  single  arm  mode,  the  two-arm  Michelson  mode,  the  three-arm 
Michelson  mode,  and  the  three-arm  Sagnac  mode.  The  calculated  averaged  response  functions 
are  presented  as  Figures  16  -  19. 


Averaged  Response  in  Single  Arm  Mode 


-  GR  Polarizations 

-  Scalar  Polarization 


Log  (u) 


Figure  16:  The  single  arm  response. 


The  comparative  response  functions  are  built  on  the  assumption  that  the  amplitude  of  the 
scalar  polarization  is  equal  to  the  amplitude  of  the  GR  polarizations,  although  most  theories 
predict  a  smaller  amplitude  for  the  scalar  wave.  Figure  16  indicates  the  averaged  scalar  response 
is  stronger  than  the  averaged  GR  response  in  the  single  arm  mode.  The  stronger  scalar  response 
is  a  direct  result  of  the  symmetry  of  the  scalar  polarization,  which  results  in  fewer  null  responses 
associated  with  the  source  direction  of  the  gravitation  wave. 
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Averaged  Response  in  Two-Arm  Michelson  Mode 
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Figure  17:  The  two-arm  Michelson  response. 


The  two-arm  Michelson  mode  enhances  the  response  to  GR  polarizations  for  low 
frequencies  while  suppressing  the  scalar  response.  The  results  are  intuitively  clear  considering 
GR  polarizations  act  to  stretch  one  arm  while  shrinking  the  second  causing  a  greater  difference 
in  arm-lengths.  The  scalar  polarization  in  low  frequencies  acts  to  stretch  or  shrink  both  arms 
simultaneously,  thereby  minimizing  the  difference  in  arm-lengths.  A  similar  argument  can  be 
posed  for  the  three-arm  Michelson  mode  (See  Fig.  18).  However,  it  is  also  important  to  note  that 
at  high-frequencies  for  which  the  wavelengths  of  gravitational  waves  are  shorter  than  the  arm 
length  of  the  detector,  the  periodic  cancellations  that  reduce  the  response  to  the  GR  polarizations 
provide  a  relative  enhancement  for  the  response  to  the  scalar  polarization.  As  the  program 
calculates  the  response  function  it  assumes  an  initial  phase  of  the  gravitational  wave  at  the 
detector.  Because  the  phase  of  the  gravitational  wave  producing  the  maximum  response  is 
unknown,  the  program  increments  this  phase,  finding  and  storing  the  maximum,  which  is  the 
amplitude  of  the  sinusoidal  response. 


(toajT) 
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Averaged  Michelson  Three  Arm  Mode 


GR  Polarizations 
Scalar  Polarization 


Frequency  (Log  (u)) 


Figure  18:  The  three-arm  Michelson  response.. 
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Averaged  Response  in  Three-Arm  Sagnac  Mode 
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Figure  19:  The  three-arm  Sagnae  response. 
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The  three-arm  Sagnac  mode  provides  the  greatest  response  to  the  sealar  polarization. 
Though  the  response  funetion  suggests  the  three-arm  Sagnac  mode  is  the  most  likely  to  detect  a 
scalar  gravitational  wave,  a  complete  sensitivity  analysis  of  LISA  requires  an  examination  of  the 
noise  in  the  three-arm  Sagnac  mode,  which  may  be  significantly  greater  than  the  noise  in  the 
Michelson  modes. 

The  LISA  response  to  gravitational  waves  from  a  particular  direction  in  the  sky  will  vary 
as  the  detector  orbits  the  sun.  In  order  to  examine  the  response  of  LISA  to  gravitational  waves  of 
various  frequencies  and  polarizations  originating  from  the  center  of  the  galaxy,  a  three 
dimensional  plot  was  generated  for  each  mode  of  data  collection.  For  each  three  dimensional 
plot,  the  abscissa  was  chosen  to  indicate  orbital  position  ranging  from  0  to  360  degrees.  The 
ordinate  indicates  the  log  of  the  angular  frequency  of  the  incident  gravitational  wave  (log  (f), 
where  f  is  measured  in  Hz);  and  the  log  of  the  magnitude  of  the  response  is  plotted  on  the  z  axis. 
The  comparative  plots  for  each  mode  of  data  collection  are  presented  as  Figures  20  -  24. 

Figure  20  displays  the  single  arm  response  for  each  arm  of  LISA  to  GR  gravitational 
waves  originating  from  the  direction  of  the  galactic  center.  The  fluctuations  in  the  response 
throughout  the  orbit  depend  heavily  upon  the  initial  orientation  of  each  arm.  Li  exhibits  a  fairly 
uniform  response  across  the  orbital,  so  it  was  chosen  for  a  comparison  with  the  single-arm  mode 
response  to  scalar  gravitational  waves  (See  Fig.  21). 


Response  (Log  (R)) 
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Figure  20:  Single  arm  response  functions  for  averaged  GR  polarizations  propagating  from  the 
galactic  center. 
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Figure  21 :  Single  arm  response  functions  for  the  scalar  and  averaged  GR  polarizations 
propagating  from  the  galactic  center. 
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Figure  22:  The  two-arm  Michelson  response  functions  for  the  scalar  and  averaged  GR 
polarizations  propagating  from  the  galactic  center. 


The  response  functions  plotted  for  the  two-arm  Michelson  mode  (See  Fig.  22)  indicate 
that  the  maximum  responses  to  low-frequency  scalar  gravitational  waves  will  occur  at  orbital 
positions  that  differ  from  those  exhibiting  the  maximum  responses  to  GR  polarizations.  The 
separation  of  maximum-response  locations  provides  a  distinguishing  feature  when  the  source 


location  is  known. 
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Figure  23:  The  three-arm  Miehelson  response  flinetions  for  the  sealar  and  averaged  GR 
polarizations  propagating  from  the  galaetie  eenter. 

The  sharp  dips  in  the  LISA  response  to  scalar  waves  measured  in  the  three-arm 
Michelson  mode  provide  an  essential  distinguishing  feature  between  the  scalar  and  GR 
polarizations  (See  Fig.  23).  A  dramatic  drop  in  signal  strength  at  four  locations  in  LISA’s  orbit 
would  indicate  the  presence  of  a  scalar  polarized  gravitational  wave.  Moreover,  the  location  of 
signal  losses  provides  information  about  the  source  direction.  Comparing  signal  response 
fluctuations  in  the  two-arm  Michelson  mode  with  those  measured  in  the  three-arm  Michelson 


mode  will  provide  greater  precision  in  the  determination  of  polarization  and  source  direction. 
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Figure  24:  The  three-arm  Sagnac  response  functions  for  the  scalar  and  averaged  GR  polarizations 
propagating  from  the  galactic  center. 


The  orbital  variations  in  response  displayed  in  the  three-arm  Sagnac  mode  are  similar  to 
those  produced  in  the  two-arm  Michelson  mode;  however,  the  scalar  response  exhibits  a 
significant  enhancement  in  magnitude.  The  Sagnac  mode  provides  improved  scalar  response 
because  the  lengths  of  the  arms  are  summed  rather  than  subtracted  from  one  another. 

To  examine  the  LISA  response  to  gravitational  waves  of  a  particular  frequency 
originating  from  the  direction  of  the  galactic  center,  data  from  the  low  frequency  end  of  the  three 
dimensional  plots  were  converted  to  polar  form.  The  magnitude  of  the  response  was  scaled  by  a 
factor  of  100  before  its  logarithm  was  computed  so  that  the  logarithm  would  be  positive  and  then 
plotted  as  the  radial  coordinate  of  a  polar  plot.  The  angular  coordinate  in  the  polar  plots 
corresponds  to  the  orbital  position  of  LISA  in  the  heliocentric  frame.  Two  polar  graphs  were 
created  to  display  the  response  functions  for  GR  polarizations  and  the  scalar  polarization  for  each 
mode  of  data  collection.  This  representation  of  the  response  as  a  function  of  orbital  position 


provides  characteristic  patterns  that  distinguish  each  polarization  and  mode  combination.  The 
plots  are  presented  as  Figures  25  and  26. 
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Averaged  Einstein  Polarization  Responses  in  Several  Data  Collection  Modes 


Figure  25:  The  response  functions  calculated  with  respect  to  the  polarizations  predicted  by  GR  for 
four  different  modes  of  data  collection. 
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Scalar  Polarization  Responses  in  Several  Data  Collection  Modes 


Figure  26:  The  response  funetions  calculated  with  respect  to  the  scalar  polarization  for  four 
different  modes  of  data  collection. 

It  is  important  to  note  that  the  response  functions  were  generated  assuming  that  the  three 
arms  of  LISA  are  of  equal  length  and  the  original  distance  between  pairs  of  spacecraft  is  very 
precisely  known.  The  study  published  by  LHH  presents  the  argument  that  assuming  the  arms  are 
of  equal  and  known  length  will  lead  to  response  functions  that  are  accurate  first  order 
representations  of  the  actual  unequal  arm  design  of  the  LISA  detector.  The  approximation  is 
appropriate  for  this  study,  as  it  is  mainly  concerned  with  comparing  the  LISA  response  functions 
generated  with  respect  to  the  GR  and  scalar  polarizations. 
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Conclusions 


The  primary  goal  of  this  study  was  to  determine  if  a  gravitational  wave  of  sealar 
polarization  with  suffieient  amplitude  eould  be  detected  and  distinguished  from  the  plus  and 
cross  polarizations.  The  results  of  the  study  indicate  that  LISA  will  be  able  to  detect  and  identify 
the  scalar  polarization  provided  the  amplitude  of  gravitational  waves  exhibiting  the  scalar 
polarization  is  comparable  to  the  amplitude  of  gravitational  waves  exhibiting  GR  polarizations. 
Additionally,  the  results  show  that  a  comparison  of  the  variations  in  responses  measured  by 
several  data  collection  modes  will  indicate  the  direction  of  the  source.  When  LISA  becomes 
operational,  variances  in  the  strength  of  a  detected  signal  may  be  compared  to  signature  features 
of  the  response  functions,  such  as  the  number  of  dips  and  their  relative  sharpness,  in  order  to 
determine  the  propagation  direction  and  polarization  of  the  measured  gravitational  wave.  These 
results  will  provide  input  for  future  researchers  to  determine  the  optimal  data  collection  method 
to  detect  gravitational  waves  of  a  given  polarization,  frequency,  and  source  direction  and  may 
also  be  used  to  develop  algorithms  to  identify  the  polarizations  and  directionality  of  gravitational 


waves  from  data  transmitted  from  the  detector. 
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Appendix  A 

A  major  task  in  developing  the  response  function  was  to  describe  the  gravitational  wave 
in  the  LISA  frame.  The  challenge  in  determining  the  description  of  the  gravitational  wave  in  the 
LISA  frame  resides  in  the  fact  that  LISA’s  orientation  relative  to  the  source  changes  as  it  orbits 
the  sun.  Several  orbital  parameters  outlined  in  the  Pre-phase  A  Report  [7]  developed  by  NASA 
made  it  possible  to  include  the  orbital  dynamics  in  the  mathematical  model  constructed  within 
this  study.  The  following  parameters  were  included  in  the  model. 

First,  the  plane  formed  by  the  three  spacecraft  forms  a  60  degree  angle  with  the  x-y  plane 
of  the  heliocentric  frame  (See  Fig.  Al).  Secondly,  the  equilateral  triangle  formed  by  the 
spacecraft  rotates  360  degrees  clockwise  about  the  z  axis  of  the  LISA  frame  as  it  follows  earth  in 
one  revolution  about  the  sun  (See  Fig.  Al). 


Figure  Al :  An  illustration  of  the  60  degree  inelination  to  the  plane  of  the  eeliptie,  and  the 
eloekwise  rotation  of  LISA  over  the  period  of  the  orbit  [10]. 

Third,  the  orbital  position  of  LISA  has  been  divided  into  days,  with  day  zero  arbitrarily  chosen 
such  that  the  x  axis  of  the  LISA  frame  parallels  the  x  axis  of  the  heliocentric  frame.  The  x  axis 
of  the  heliocentric  frame  points  from  the  sun  to  the  First  Point  of  Aries. 
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Five  rotations  are  used  to  transform  from  the  deseriptions  of  the  gravitational  wave  in  the 
source  frame  to  the  description  of  the  gravitational  wave  in  the  LISA  frame.  Each  rotation  is 
described  in  terms  of  the  polar  angles  0  and  (j)  defined  in  the  heliocentric  frame  (See  Fig.  A2). 
The  first  two  rotations  bring  the  source  frame  into  the  heliocentric  frame,  while  the  last  three 
rotations  transform  between  the  heliocentric  frame  and  the  LISA  frame. 

The  y  axis  of  the  source  frame  is  assumed  to  be  coplanar  with  the  z  axis  of  the 
heliocentric  frame.  (Note:  Any  adverse  implications  of  this  assumption  are  nullified  by 
averaging  the  response  over  the  plus  and  cross  polarizations;  and  the  assumption  has  no  effect  on 
the  calculation  of  the  response  with  respect  to  the  scalar  polarization  due  to  the  symmetric 
properties  of  the  strain  pattern.)  The  rotations  are  listed  in  order  in  Table  A1 . 
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Rotation  Axis  of  Rotation  Angle  Subtended 


Ri 

X 

7t:-0 

R2 

z 

7t:/2-(j) 

R3 

z 

27t:t/365 

R4 

X 

n/3 

R5 

z 

-27t:t/365 

Table  A1 :  The  first  two  rotations  listed  are  used  to  convert  the  source  frame  into  the  heliocentric 
frame.  The  remaining  three  rotations  convert  to  the  LISA  frame.  Here  the  variable  t  represents 
the  number  of  days  elapsed  in  LISA’s  orbit  about  the  sun. 


Rotations  about  the  x  and  z  axes  are  provided  by  the  following  equations. 


R 


X 


1  0  0 
0  cos(^ir)  sin(^) 

0  -  sin(^)  oos(^) 


R 


z 


oos(^)  sin(^)  0 
-  sin(^)  cos(^)  0 
0  0  1 


(Al) 


(A2) 


Here  \\>  is  the  angle  subtended  in  the  rotation.  By  eonvention,  the  wave  propagates  in  the  Z 
direetion  of  the  souree  frame.  So  the  propagation  direetion  of  the  gravitational  wave  in  the  LISA 
frame  (klisa),  is  given  by  equation  A3. 


^LISA 


0 

0 

1 


(A3) 
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The  metric  describing  the  gravitational  wave  represents  a  tensor,  so  the  original  matrix  must  also 
be  multiplied  by  the  transpose  of  each  rotation. 

WXR^^R'M  (A4) 


The  perturbation  matrix  in  the  LISA  frame  is  therefore  written  as: 


K^,,,={R.RAR2RlH,,_]RlR‘AKR^)^^os{e+d)  (as) 

These  rotations  are  construction  on  the  assumption  that  the  initial  orbital  position  of  LISA  is  on 
the  y-axis  of  the  heliocentric  frame. 


42 


Appendix  B 

The  following  MATLAB®  code  calculates  the  response  functions  averaged  over  all 

possible  source  directions. 


function  Z  =  AverageAllSky  () 

%  Function  that  will  determine  the  LISA  response  for  plus,  cross,  and 
%  scalar  polarizations  in  an  all  sky  average.  The  program  also  allows  for  a 
%  linear  combination  of  cross  and  plus,  and  a  linear  combination  of  cross 
%  and  scalar  polarizations. 

%  Begin  clock  to  measure  the  program  run-time, 
tic 


%  The  perturbation  matricies 
hcross  =  [0  1  0;  1  0  0;  0  0  0]; 
hplus  =[10  0;  0  -1  0;  0  0  0]; 
hscalar  =  [1  0  0;  0  1  0;  000]; 

%  The  Mixfactor  is  the  coefficient  in  a  linear  combination  of  cross  and 
%  scalar  polarizations.  For  example:  Mixfactor  =  .1  would  mean  that  the 
%  matrix  element  of  the  calculations  consists  of  hx  +  .l*hs,  where  hx  is 
%  the  matrix  element  for  the  cross  polarization  and  hs  is  the  matrix 
%  element  for  the  scalar  polarization. 

Mixfactor  =  .1; 

%  tstep  gives  the  sample  size  in  days.  For  example:  tstep  =1  then  the 
%  response  will  be  calculated  for  each  day  of  its  orbit, 
tstep  =  10; 

tdays  =  [ 0 : tstep : 365] ; 

%  The  range  of  frequencies 
uu  =  2*pi* (10 . ^  (-4 :  .  01 : 1 . 2)  )  ; 

%  Matricies  used  to  store  values 
[r,  c]  =  size  (uu)  ; 

[rr,cc]  =  size  (tdays); 

A  =  ones (c, cc)  ; 

B  =  A; 

C  =  A; 

D  =  A; 

E  =  A; 

F  =  A; 


ScalarA 

=  A 

ScalarB 

=  A 

ScalarC 

=  A 

ScalarD 

=  A 

ScalarE 

=  A 

ScalarF 

=  A 

Combo A  = 

A; 

ComboB  = 

A; 

ComboC  = 

A; 

ComboD  = 

A; 

ComboE  = 

A; 

ComboF  = 

A; 

%  The  values  associated  with  alpha  in  the  LISA  frame  and  store  in 
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%  memory  for  later  use.  Alpha  is  the  angle  between  an  arm  of  LISA  and  the 

%  x-axis  of  the  LISA  frame. 

csal=cos (pi/12) ; 

csalS=csal^2 ; 

csa2=cos (5*pi/12)  ; 

csa2S=csa2^2 ; 

csa3=cos ( 3*pi/4 )  ; 

csa3S=csa3^2 ; 

snal  =  sin (pi/12 )  ; 

snalS=snal''2  ; 

sna2  =  sin (5*pi/12)  ; 

sna2S=sna2^2 ; 

sna3=sin (3*pi/4 )  ; 

sna3S=sna3^2 ; 

%  Counting  numbers  used  to  designate  where  data  is  stored. 
n=0; 

N=0; 

nn=0; 

%  The  ranges  and  step  sizes  for  theta,  phi,  and  epsilon, 
thetal  =  [pi/16  :  (pi/8 ):  (ll*pi/l 6) ] ; 

[rtheta,  ctheta]  =  size  (thetal); 

phil  =  [pi/16:  (pi/8)  :  (31*pi/16) ] ; 

[rphi,  cphi]  =  size  (phil); 

epsilonstep  =  1/16; 

epsilonl  =  [  (  (pi*epsilonstep) 12)  :  (pi*epsilonstep)  :  (  (2/epsilonstep) -1) *pi. . . 

* (epsilonstep/2 ) ] ; 

[reps, ceps]  =  size  (epsilonl); 

%  The  final  rotation  that  orients  the  LISA  frame  sixty  degrees 
%  relative  to  the  plane  of  the  ecliptic. 

RorbitTilt  =[10  0;  0  cos  (pi/3)  sin  (pi/3);  0  -sin  (pi/3)  cos  (pi/3)]; 


%  A  series  of  nested  loops  numerically  evaluates  the  response  function, 
for  theta  =  thetal 

%  The  values  associated  with  theta  including  the  first  rotation 
%  required  to  translate  between  the  source  frame  and  the  heliocentric  frame. 
cth=  cos (theta); 
sth=  sin (theta); 

RTheta  =[10  0;  0  cos (pi-theta)  sin (pi-theta) ;  0  -sin (pi-theta) .. . 

cos  (pi-theta) ] ; 

for  phi  =  phil 

%  Each  value  of  nn  corresponds  to  a  unique  source  location. 
nn=  nn+1; 

%  The  values  used  in  the  calculation  of  the 
%  response  function  that  include  the  variable  phi. 
muone=  -sth*cos (phi-pi/ 12 ) ; 
mutwo=  -sth*cos (phi-5*pi/12 )  ; 
muthree=  -sth*cos (phi-3*pi/ 4 )  ; 

%  The  second  rotation  required  to  translate  between  the 
%  source  and  heliocentric  frame. 

RPhi  =  [cos (pi/2-phi)  sin (pi/2-phi)  0;  -sin (pi/2-phi)  cos (pi/2-phi) . . . 

0 ;  0  0  1  ]  ; 
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%  The  perturbation  matrices  in  the  heliocentric  frame. 
hHcross  =  RPhi ' *RTheta ' *hcross*RTheta*RPhi; 
hHplus  =  RPhi ' *RTheta ' *hplus*RTheta*RPhi ; 
hHscalar  =  RPhi ' *RTheta ' *hscalar*RTheta*RPhi ; 

%  Clear  the  matrices  used  to  temporarily  store  the  response  function 
%  calculated  for  a  particular  value  of  theta  and  phi. 

A  =  zeros (c, cc) ; 

B  =  A; 

C  =  A; 

D  =  A; 

E  =  A; 

F  =  A; 

ScalarA  =  A; 

ScalarB  =  A; 

ScalarC  =  A; 

ScalarD  =  A; 

ScalarE  =  A; 

ScalarF  =  A; 

ComboA  =  A; 

ComboB  =  A; 

ComboC  =  A; 

ComboD  =  A; 

ComboE  =  A; 

ComboF  =  A; 

%  Reset  counting  variable  used  to  indicate  frequency. 

N=0; 

for  u  =  uu 

%  N  and  n  are  counting  numbers  that  control  the  storage  of 
%  values  in  Matrices  A,  B,  and  C.  N  corresponds  to  the 
%  frequency  and  n  corresponds  to  the  orbital  position. 

N  =  N  +  1; 
n=0; 

for  t  =  tdays 
n=n+l ; 

%  Rotation  matrices  that  account  for  EISA's  orbit  about 
%  the  sun. 

RphiLISA  =  [cos  (  ( t*2*pi ) / 365 )  sin  (  ( t*2*pi ) /365 )  0;... 

-sin  (  (t*2*pi) 7365)  cos  (  ( t*2 *pi ) 73 65 )  0;  0  0  1]; 
RbackRotate  =  [cos  (- (t*2*pi) 7365)  sin  (- (t*2*pi) 7365)  0;... 

-sin  (- (t*2*pi) 7365)  cos  (- ( t*2 *pi ) 73 65 )  0;  0  0  1]; 

%  The  perturbation  matrices  described  in  the 
%  LISA  frame.  Three  Euler  angles  are  used  to  transform 
%  the  perturbation  matrices  described  in  the 
%  heliocentric  frame  to  those  described  in  the  LISA 
%  frame. 

hLcross  =  RbackRotate*RorbitTilt*RphiLISA*hHcross*RphiLISA' . . . 
*RorbitTilt ' *RbackRotate ' ; 

hLplus  =  RbackRotate*RorbitTilt*RphiLISA*hHplus*RphiLISA ' . . . 
*RorbitTilt ' *RbackRotate ' ; 

hLscalar=  RbackRotate*RorbitTilt*RphiLISA*hHscalar*RphiLISA ' . . . 
*RorbitTilt ' *RbackRotate ' ; 


%  The  time  independent  components  associated  with  each 
%  response  function  evaluated. 

CrossMetl=  . 5* ( (hLcross ( 1 , 1 ) *csalS )  +  (hLcross ( 1 , 2 ) *csal* snal ) . 

+  (hLcross (2, 1) *snal*csal)  +  (hLcross (2 , 2 ) * snalS ) )  ; 
CrossMet2=  . 5* (  (hLcross  ( 1 , 1 ) *csa2S )  +  (hLcross ( 1 , 2 ) *csa2* sna2 )  . 

+  (hLcross (2, 1) *sna2*csa2)  +  (hLcross  (2 , 2 ) * sna2S ) ) ; 
CrossMet3=  . 5* (  (hLcross ( 1 , 1 ) *csa3S )  +  (hLcross ( 1 , 2 ) *csa3* sna3 )  . 
+  (hLcross (2, 1) *sna3*csa3)  +  (hLcross (2 , 2 ) * sna3S ) ) ; 

PlusMetl=  . 5*  (  (hLplus  (1, 1) *csalS)  +  (hLplus (1 , 2 ) *csal*snal )  . . . 

+  (hLplus  (2, 1) *snal*csal)  +  (hLplus (2 , 2 ) *snalS) ) ; 

PlusMet2=  . 5* ( (hLplus (1, 1) *csa2S)  +  (hLplus  (1 , 2 ) *csa2*sna2 ).. . 

+  (hLplus (2, 1) *sna2*csa2)  +  (hLplus (2 , 2 ) * sna2S ) ) ; 

PlusMet3=  . 5* ( (hLplus (1, 1) *csa3S)  +  (hLplus (1 , 2 ) *csa3*sna3) .. . 

+  (hLplus (2, 1) *sna3*csa3)  +  (hLplus (2 , 2 ) * sna3S ) ) ; 

ScalarMetl=  . 5* ( (hLscalar  ( 1 , 1 ) *csalS )  +  (hLscalar  ( 1 , 2 ) *csal .  .  . 

*snal)  +  (hLscalar (2, 1) *snal*csal)  +  (hLscalar (2, 2) *snalS) ) 
ScalarMet2=  . 5*  (  (hLscalar ( 1 , 1 ) *csa2S )  +  (hLscalar ( 1 , 2 ).. . 

*csa2*sna2)  +  (hLscalar (2 , 1 ) *sna2*csa2 )  + (hLscalar  (2 , 2 ).. . 
*sna2S) ) ; 

ScalarMet3=  . 5*  (  (hLscalar ( 1 , 1 ) *csa3S )  +  (hLscalar ( 1 , 2 ) *csa3 .. . 
*sna3)  +  (hLscalar (2, 1) *sna3*csa3)  +  (hLscalar ( 2 , 2 ) *sna3S ) ) 

%  Resetting  variables  used  in  calculating  the  response 
%  function  to  zero  for  the  next  iteration. 

DXl sum=0 ; 

DX2  sum=0 ; 

DX3sum=0; 

DX12sum=0 ; 

DX123Msum=0; 

DX123Ssum=0; 

DXScalarl sum=0 ; 

DXScalar2  sum=0 ; 

DXScalar3sum=0; 

DXScalarl2sum=0  ; 

DXScalarM12  3sum=0 ; 

DXScalarS12  3sum=0 ; 

DXC  omb  o 1 s  um= 0 ; 

DXC  omb  o  2  s  um= 0 ; 

DXCombo3sum=0 ; 

DXCombol2  sum=0 ; 

DXComboMl 2  3  sum=0 ; 

DXComboS123sum=0; 

%  Averaging  over  epsilon  accounts  for  rotation  about 
%  the  third  Euler  angle  required  in  transforming  from 
%  the  source  frame  to  the  heliocentric  frame, 
for  epsilon  =  epsilonl 

%  This  combination  of  the  plus  and  cross  metric 
%  components  accounts  for  the  third  rotation 
%  mentioned  in  the  comment  on  line  211  (directly 
%  above  epsilon  for  loop) . 

metepsone  =  cos (2*epsilon) *CrossMetl  -  sin (2*epsilon) . . . 
*PlusMetl ; 

metepstwo  =  cos  (2*epsilon) *CrossMet2  -  sin (2*epsilon)  . . . 
*PlusMet2 ; 

metepsthree  =  cos (2*epsilon) *CrossMet3  - 
*PlusMet3 ; 


sin (2*epsilon) . . . 


%  Combines  the  scalar  polarization  with  the  mixture 
%  of  GR  plus  and  cross  polarizations. 

ComboCrossScalarl  =  metepsone  +  Mixf actor *ScalarMetl ; 
ComboCrossScalar2  =  metepstwo  +  Mixf actor *ScalarMet2 ; 
ComboCrossScalarS  =  metepsthree  +  Mixf actor *ScalarMet3 ; 

%  Resetting  values  used  in  delta  loop. 

DX1=0; 

DX2=0; 

DX3=0; 

DX12=0; 

DX123M=0; 

DX123S=0; 

DXScalarl=0; 

DXScalar2=0; 

DXScalar3=0 ; 

DXScalarl2=0; 

DXScalarM123=0; 

DXScalarS123=0; 

DXCombol=0; 

DXCombo2=0; 

DXCombo3=0 ; 

DXCombol2=0; 

DXComboM123=0; 

DXComboS123=0; 

%  Delta  is  the  phase  constant  for  the  incident 
%  wave.  ( kx-wt+delta)  This  portion  of  the  program 
%  finds  the  amplitude  of  the  response.  The 
%  response  will  vary  depending  upon  the  initial 
%  phase  of  the  gravitational  wave.  As  the  program 
%  cycles  through  possible  values  for  the  initial 
%  phase,  it  stores  the  maximum  response, 
for  delta  =  pi/ 32 : pi/ 1 6 : 31*pi/32 ; 

%  fl,  f2,  and  f3  include  the  time  dependent  parts  for 
%  the  response  function. 

fl=  ( (sin (delta-muone*u) -sin (delta-u) ) / (u* (1-muone) ) . 

+  (sin  (u+delta) -sin (delta-muone*u) ) / (u* (1+muone) ) ) 
f2=  (  (sin (delta-mutwo*u) -sin (delta-u) ) / (u* (1-mutwo) )  . 

+ (sin (u+delta) -sin (delta-mutwo*u) ) / (u* (1+mutwo) ) ) 
f3=  ( (sin (u* (1+muone) -delta) -sin (mutwo*u-delta) ) . . . 

/  (u* (1-muthree) )  +  (sin (u* (1-muone) +delta)  . . . 

+sin (mutwo*u-delta) ) / (u* ( 1+mu three) ) ) ; 

%  deltaX  values  correspond  to  the  change  in  length  of 
%  each  arm.  This  combines  the  time  independent 
%  components  with  the  time  dependent  component  to 
%  achieve  the  complete  response  of  LISA  to  the 
%  gravitational  wave.  XI  corresponds  to  LI,  X2  to  L2, 
%  and  so  on.  Similar  calculations  were  computed  for 
%  the  scalar  response  and  the  response  to  the 
%  combination  of  scalar  and  GR  polarizations. 
deltaXl=  metepsone* fl ; 
deltaX2=  metepstwo* f 2 ; 
deltaX3=  metepsthree*f3; 

%  In  order  to  calculate  the  response  we  must  find  the 
%  value  of  delta  that  yields  the  maximum  result 
%  (indicating  amplitude  of  response)  . 
if  abs(deltaXl)  >  DXl 


DXl  =  abs (deltaXl ) ; 


end 

if  abs(deltaX2)  >  DX2 
DX2  =  abs (deltaX2 ) ; 

end 

if  abs(deltaX3)  >  DX3 
DX3  =  abs (deltaX3 )  ; 

end 

if  abs (deltaXl-deltaX2 )  >  DX12 
DX12  =  abs (deltaXl-deltaX2 ) ; 

end 

if  abs (deltaXl+deltaX2-2*deltaX3)  >  DX123M 
DX123M  =  abs (deltaXl+deltaX2-2*deltaX3) ; 

end 

if  abs (deltaXl+deltaX2+deltaX3)  >  DX123S 
DX123S  =  abs (deltaXl+deltaX2+deltaX3) ; 

end 

deltaXScalarl  =  ScalarMetl*fl ; 
deltaXScalar2  =  ScalarMet2*f2 ; 
deltaXScalar3  =  ScalarMet3*f3; 

if  abs (deltaXScalarl )  >  DXScalarl 

DXScalarl  =  abs (deltaXScalarl ) ; 

end 

if  abs (deltaXScalar2)  >  DXScalar2 
DXScalar2  =  abs (deltaXScalar2 ) ; 

end 

if  abs (deltaXScalarl)  >  DXScalarl 
DXScalarl  =  abs (deltaXScalarl ) ; 

end 

if  abs (deltaXScalarl-deltaXScalar2 ) >  DXScalarl! 

DXScalarl!  =  abs (deltaXScalarl-deltaXScalar! ) ; 

end 

if  abs  (deltaXScalarl+deltaXScalar2-2 *deltaXScalar3 ) 

>  DXScalarM123 

DXScalarM123  =  abs  (deltaXScalarl+deltaXScalar! . 
-2*deltaXScalar3) ; 

end 

if  abs  (deltaXScalarl+deltaXScalar2+deltaXScalar3 ) . . 

>  DXScalarSl!! 

DXScalarS123  =  abs  (deltaXScalarl+deltaXScalar!. 
tdeltaXScalarl) ; 

end 

deltaXCombol  =  ComboCrossScalarl*fl ; 
deltaXCombo!  =  ComboCrossScalar!*f! ; 
deltaXCombol  =  ComboCrossScalarl* f 3 ; 

if  abs (deltaXCombol )  >  DXCombol 
DXCombol  =  abs (deltaXCombol ) ; 

end 

if  abs (deltaXCombo! )  >  DXCombo! 

DXCombo!  =  abs (deltaXCombo! ) ; 

end 

if  abs (deltaXCombol )  >  DXCombol 
DXCombol  =  abs (deltaXCombol ) ; 

end 

if  abs (deltaXCombol-deltaXCombo! ) >  DXCombol! 

DXCombol!  =  abs (deltaXCombol-deltaXCombo! ) ; 

end 

if  abs  (deltaXCombol+deltaXCombo!-!*deltaXCombo3) . . . 


>  DXComboM123 

DXComboM123  =  abs  (deltaXCombol+deltaXCombo2 . 
-2*deltaXCombo3 ) ; 

end 

if  abs  (deltaXCombol+deltaXCombo2+deltaXCombo3) . . 

>  DXComboS123 

DXComboS123  =  abs  (deltaXCombol+deltaXCombo2+ 
deltaXCombo3) ; 

end 


end  %  End  delta  loop. 

%  Add  the  response  functions  for  each  arm  and 
%  each  combination  of  arms. 

DX12sum  =  DX12sum  +  abs(DX12); 

DXlsum  =  DXlsum  +  abs(DXl); 

DX2sum  =  DX2sum  +  abs(DX2); 

DX3sum  =  DX3sum  +  abs(DX3); 

DX123Msum  =  DX123Msum  +  abs (DX123M) ; 

DX123Ssum  =  DX123Ssum  +  abs (DX123S) ; 

DXScalarlsum=DXScalarlsum  +  abs (DXScalarl ) ; 
DXScalar2sum=DXScalar2sum  +  abs (DXScalar2 ) ; 
DXScalar3sum=DXScalar3sum  +  abs (DXScalar3 ) ; 
DXScalarl2sum=DXScalarl2sum  +  DXScalarl2; 
DXScalarM123sum=DXScalarM123sum  +  DXScalarM12 3 ; 
DXScalarS123sum=DXScalarS123sum  +  DXScalarS123; 


DXCombolsum=DXCombolsum  +  abs (DXCombol ) ; 
DXCombo2sum=DXCombo2sum  +  abs  (DXCombo2 ) ; 
DXCombo3sum=DXCombo3sum  +  abs (DXCombo3) ; 
DXCombol2sum=DXCombol2sum  +  DXCombol2; 
DXComboM123sum=DXComboM123sum  +  DXComboM123 ; 
DXComboS123sum=DXComboS123sum  +  DXComboS123 ; 


end  %  End  epsilon  loop 

%  Stores  the  sum  of  the  responses  for  each  angle  sampled. 
A(N,n)  =  DXlsum; 

B  (N, n)  =  DX2sum; 

C  (N, n )  =  DX3sum; 

D(N,n)  =  DX12sum; 

E(N,n)  =  DX123Msum; 

F(N,n)  =  DX123Ssum; 


ScalarA (N, n) 
ScalarB  (N, n) 
ScalarC (N, n) 
ScalarD  (N, n) 
ScalarE  (N, n) 
ScalarF  (N, n) 


DXScalarlsum; 

DXScalar2sum; 

DXScalarlsum; 

DXScalarl2sum; 

DXScalarM123sum; 

DXScalarS123sum; 


ComboA (N, n) 
ComboB (N, n) 
ComboC (N, n) 
ComboD (N, n) 
ComboE (N, n) 
ComboF (N, n) 


DXCombolsum; 
DXCombo2sum; 
DXCombolsum; 
DXCombol2  sum; 
DXComboMl 2  3  sum; 
DXComboS123sum; 


end  %end  t  loop 
end  %end  u  loop 


In  order  to  achieve  an  equally  weighted  average  over  all  possible 
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%  source  directions  some  weighting  factors  had  to  be  included.  The 
%  sin  (theta)  factor  accounts  for  the  fact  that  on  a  (lat,  long) 

%  grid,  such  as  the  (theta,  phi)  coordinates  chosen,  there  will  be 
%  more  functions  sampled  near  the  poles  than  near  the  equatorial 
%  plane. 

aA(;,;,nn)  =  A*sth; 
aB ( ; , : , nn)  =  B*sth; 
aC ( : , : , nn)  =  C*sth; 
aD ( ; , ; , nn)  =  D*sth; 
aE ( ; , : , nn)  =  E*sth; 
aF ( ; , : , nn)  =  F*sth; 

aScalarA ( ; , ; , nn )  =  ScalarA*sth; 
aScalarB ( ; , ; , nn)  =  ScalarB*sth; 
aScalarC  (  ; ,  : , nn)  =  ScalarC*sth; 
aScalarD ( ; , ; , nn)  =  ScalarD*sth; 
aScalarE  (  ; ,  ; , nn)  =  ScalarE*sth; 
aScalarF ( ; , : , nn)  =  ScalarF*sth; 

aComboA ( ; , ; , nn)  =  ComboA*sth; 
aComboB ( ; , ; , nn)  =  ComboB*sth; 
aComboC ( : , : , nn)  =  ComboC*sth; 
aComboD ( ; , : , nn)  =  ComboD*sth; 
aComboE ( ; , ; , nn)  =  ComboE*sth; 
aComboF ( ; , ; , nn)  =  ComboF*sth; 

end  %  End  phi  loop 

end  %  End  theta  loop 

%  Sum  the  response  functions  and  divide  by  the  respective  number  of 
%  components  summed.  Note  the  three  dimensional  matrix  must  first  be  summed 
%  along  one  axis,  then  the  other. 

RTaA  =  sum ( aA, 3 ) / ( 4 *pi*ctheta*cphi*ceps ) ; 

RTaB  =  sum ( aB, 3 ) / ( 4 *pi*ctheta*cphi*ceps ) ; 

RTaC  =  sum(aC, 3) / (4*pi*ctheta*cphi*ceps)  ; 

RTaD  =  sum ( aD, 3 ) / ( 4 *pi*ctheta*cphi*ceps ) ; 

RTaE  =  sum ( aE, 3 ) / ( 4 *pi*ctheta*cphi*ceps ) ; 

RTaF  =  sum ( aF, 3 ) / ( 4 *pi*ctheta*cphi*ceps ) ; 

RTaScalarA  =  sum (aScalarA, 3 )/( 4 *pi*ctheta*cphi*ceps )  ; 

RTaScalarB  =  sum (aScalarB, 3 )/( 4 *pi*ctheta*cphi*ceps ) ; 

RTaScalarC  =  sum ( aScalarC, 3 ) / ( 4 *pi*ctheta*cphi*ceps ) ; 

RTaScalarD  =  sum ( aScalarD, 3 )/( 4 *pi*ctheta*cphi*ceps )  ; 

RTaScalarE  =  sum (aScalarE, 3 )/( 4 *pi*ctheta*cphi*ceps ) ; 

RTaScalarF  =  sum (aScalarF, 3 )/( 4 *pi*ctheta*cphi*ceps ) ; 

RTaComboA  =  sum (aComboA, 3 ) / ( 4 *pi*ctheta*cphi*ceps ) ; 

RTaComboB  =  sum (aComboB, 3 ) / ( 4 *pi*ctheta*cphi*ceps ) ; 

RTaComboC  =  sum (aComboC, 3 ) / ( 4 *pi*ctheta*cphi*ceps ) ; 

RTaComboD  =  sum (aComboD, 3 ) / ( 4 *pi*ctheta*cphi*ceps ) ; 

RTaComboE  =  sum (aComboE, 3 ) / ( 4 *pi*ctheta*cphi*ceps ) ; 

RTaComboF  =  sum (aComboF, 3 ) / ( 4 *pi*ctheta*cphi*ceps ) ; 

%  Sums  along  the  second  axis,  and  divides  by  the  number  of  elements  summed. 
RTaA  =  sum (RTaA, 2 ) /cc; 

RTaB  =  sum (RTaB, 2 ) /cc; 

RTaC  =  sum (RTaC, 2 ) /cc; 

RTaD  =  sum (RTaD, 2 ) /cc; 

RTaE  =  sum (RTaE, 2 ) /cc; 

RTaF  =  sum (RTaF, 2 ) /cc; 

RTaScalarA  =  sum (RTaScalarA, 2 ) /cc; 
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RTaScalarB  =  sum (RTaScalarB, 2 ) /cc; 

RTaScalarC  =  sum (RTaScalarC, 2 ) /cc; 

RTaScalarB  =  sum (RTaScalarB, 2 ) /cc; 

RTaScalarB  =  sum (RTaScalarB, 2 ) /cc; 

RTaScalarB  =  sum (RTaScalarB, 2 ) /cc; 

RTaComboA  =  sum (RTaComboA, 2 ) /cc; 

RTaComboB  =  sum (RTaComboB, 2 ) /cc; 

RTaComboC  =  sum (RTaComboC, 2 ) /cc; 

RTaComboB  =  sum (RTaComboB, 2 ) /cc; 

RTaComboB  =  sum (RTaComboB, 2 ) /cc ; 

RTaComboB  =  sum (RTaComboB, 2 ) /cc; 

%  Bnd  program  run  time  and  display  in  the  command  prompt, 
toe 

%  Plot  one  figure  to  indicate  the  program  has  finished  running  and  to 
%  provide  an  initial  plot  to  verify  results, 
a  =  logic (uu) ; 
figure 

plot (a,  logic ( (RTaB' ) /2) ) 

title  ('Averaged  Michelson  Three  Arm  Response') 

%  Save  the  variable  space  so  additional  plots  can  be  generated  with  the  data 
%  without  rerunning  the  program, 
save  ( ' AllSkyCCMAY ' ) ; 
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function  Z  =  ThreeDplots  () 

%  Function  that  will  determine  the  LISA  response  for  plus,  cross,  and 
%  scalar  polarizations  originating  from  the  galactic  center.  The  resulting 
%  three  dimensional  surface  plots  illustrate  the  instrument  response  with 
%  respect  to  orbital  position  and  frequency.  The  program  allows  for  a 
%  linear  combination  of  cross  and  plus,  and  a  linear  combination  of  cross 
%  and  scalar  polarizations. 

%  Begin  clock  to  measure  the  program  run-time, 
tic 

%  The  perturbation  matricies. 
hcross  =  [0  1  0;  1  0  0;  0  0  0 ] ; 
hplus  =[10  0;  0  -1  0;  0  0  0]; 
hscalar  =  [1  0  0;  0  1  0;  000]; 

%  The  Mixfactor  is  the  coefficient  in  a  linear  combination  of  cross  and 
%  scalar  polarizations.  For  example:  Mixfactor  =  .1  would  mean  that  the 
%  matrix  element  of  the  calculations  consists  of  hx  +  .l*hs,  where  hx  is 
%  the  matrix  element  for  the  cross  polarization  and  hs  is  the  matrix 
%  element  for  the  scalar  polarization. 

Mixfactor  =  .1; 

%  tstep  gives  the  sample  size  in  days.  For  example:  tstep  =1  then  the 
%  response  will  be  calculated  for  each  day  of  its  orbit, 
tstep  =  2 ; 

tdays  =  [ 0 : tstep : 365] ; 

%  Establishes  the  range  of  frequencies 
uu  =  2*pi* (10.^  (-4:  .01:  .4)  )  ; 

%  Generates  matricies  to  store  values 
[r,  c]  =  size  (uu)  ; 

[rr,cc]  =  size  (tdays); 

A  =  ones (c, cc) ; 

B  =  A; 

C  =  A; 

D  =  A; 

E  =  A; 

F  =  A; 

ScalarA  =  A; 

ScalarB  =  A; 

ScalarC  =  A; 

ScalarD  =  A; 

ScalarE  =  A; 

ScalarF  =  A; 

ComboA  =  A; 

ComboB  =  A; 

ComboC  =  A; 

ComboD  =  A; 

ComboE  =  A; 

ComboF  =  A; 


%  Counting  numbers  used  to  designate  where  data  is  stored. 
n=0; 

N=0; 

%  The  values  associated  with  alpha  in  the  LISA  frame  and  store  in 

%  memory  for  later  use.  Alpha  is  the  angle  between  an  arm  of  LISA  and  the 

%  x-axis  of  the  LISA  frame. 

csal=cos (pi/12) ; 

csalS=csal^2 ; 

csa2=cos (5*pi/12) ; 

csa2S=csa2^2 ; 

csa3=cos ( 3*pi/4 ) ; 

csa3S=csa3^2 ; 

snal=sin (pi/12 ) ; 

snalS=snal''2  ; 

sna2  =  sin ( 5*pi/12 )  ; 

sna2S=sna2^2 ; 

sna3=sin ( 3*pi/4 )  ; 

sna3S=sna3^2 ; 

%  The  ranges  and  step  size  for  epsilon, 
epsilonstep  =  1/32; 

epsilonl  =  [  (  (pi*epsilonstep) 12)  :  (pi*epsilonstep)  :  (  (2/epsilonstep) -1) *pi* . . . 

(epsilonstep/2 ) ] ; 

[reps, ceps]  =  size  (epsilonl); 

%  Rotation  matrix  for  theta 

RorbitTilt  =[10  0;  0  cos  (pi/3)  sin  (pi/3);  0  -sin  (pi/3)  cos  (pi/3) ] ; 

%  Value  for  theta  indicating  the  direction  of  the  galactic  center, 
theta  =  1.6675; 

%  Calculates  sin  and  cosine  once  and  stores  the  values  to  optimize 
%  the  program. 
cth=  cos (theta); 
sth=  sin (theta); 

%  The  rotation  matrix  for  theta.  It  rotates  (pi-theta)  about  the  x 
%  axis. 

RTheta  =[10  0;  0  cos (pi-theta)  sin (pi-theta) ;  0  -sin (pi-theta)  cos (pi-theta) ] 

%  Value  for  phi  indicating  the  direction  of  the  galactic  center, 
phi  =  4.6572; 

%  Optimizing  the  calculations  requiring  spatial  angles  theta  and  phi. 

%  Mu  values  are  necessary  in  calculating  the  time  dependant 

%  component  of  the  the  response. 

muone=  -sth*cos (phi-pi/ 12 ) ; 

mutwo=  -sth*cos  (phi-5*pi/12 ) ; 

muthree=  -sth*cos (phi-3*pi/ 4 ) ; 

%  The  rotation  matrix  for  phi.  It  rotates  pi/2-phi  about  the 
%  z  axis. 

RPhi  =  [cos  (pi/2-phi)  sin  (pi/2-phi)  0;  -sin  (pi/2-phi)  cos (pi/2-phi)  0;  0  0  1]; 

%  Generating  the  perturbation  matrices  in  the  heliocentric  frame  through  the 

%  rotations  about  theta  and  phi. 

hHcross  =  RPhi ' *RTheta ' *hcross*RTheta*RPhi; 

hHplus  =  RPhi ' *RTheta ' *hplus*RTheta*RPhi ; 

hHscalar  =  RPhi ' *RTheta ' *hscalar*RTheta*RPhi ; 

%  The  loop  used  to  calculate  the  response  and  average  over  all  sky.  u  =  uu 
%  cycles  the  calculations  for  a  range  of  frequencies 
for  u  =  uu 


o\o  o\o  D  o\o  o\o  S  o\o  o\o 


N  is  the  counting  number  that  controls  the  storage  of  response 
function  values  for  each  new  frequency. 

=  N  +  1; 

n  is  the  counting  number  that  controls  the  storage  of  response 
function  values  for  each  new  orbital  position. 

=0; 

The  t  loop  calculates  ans  stores  the  response  of  the  detector  for  each 
point  in  its  orbit  about  the  sun. 
for  t  =  tdays 

%  Increment  the  value  of  n. 
n=n+l ; 

%  The  rotation  matrices  that  depend  upon  LISA's  orbital  position. 

RphiLISA  =  [cos ( (t*2*pi) 7365)  sin ( ( t*2*pi ) / 365 )  0;  -sin (( t*2*pi ) 7365 ).. . 
cos ( (t*2*pi) 7365)  0;  001]; 

RbackRotate  =  [cos (- (t*2*pi) 7365)  sin  (- (t*2*pi) 7365)  0;  -sin  ( -  ( t*2*pi ) 7 365 ) 
cos (- (t*2*pi) 7365)  0;  0  0  1]; 

%  The  perturbation  matrices  calculated  in  the  LISA  frame. 

hLcross  =  RbackRotate*RorbitTilt*RphiLISA*hHcross*RphiLISA' *RorbitTilt ' * . . . 
RbackRotate '  ; 

hLplus  =  RbackRotate*RorbitTilt*RphiLISA*hHplus*RphiLISA '*RorbitTilt'*. . . 
RbackRotate '  ; 

hLscalar=  RbackRotate*RorbitTilt*RphiLISA*hHscalar*RphiLISA '*RorbitTilt'*. . 
RbackRotate '  ; 

%  The  time  independent  components  associated  with  each 
%  response  function  evaluated. 

CrossMetl=  . 5* ( (hLcross ( 1 , 1 ) *csalS )  +  (hLcross ( 1 , 2 ) *csal* snal )  +  ... 

(hLcross (2, 1) *snal*csal)  +  (hLcross (2 , 2 ) * snalS ) ) ; 

CrossMet2=  . 5* ( (hLcross ( 1 , 1 ) *csa2S )  +  (hLcross ( 1 , 2 ) *csa2* sna2 )  +  ... 

(hLcross (2, 1) *sna2*csa2)  +  (hLcross (2 , 2 ) * sna2S ) ) ; 

CrossMet3=  . 5* ( (hLcross ( 1 , 1 ) *csa3S )  +  (hLcross ( 1 , 2 ) *csa3* sna3 )  +  ... 
(hLcross (2, 1) *sna3*csa3)  +  (hLcross (2 , 2 ) * sna3S ) ) ; 

PlusMetl=  . 5* ( (hLplus (1, 1) *csalS)  +  (hLplus ( 1 , 2 ) *csal *snal )  +  ... 

(hLplus (2, 1) *snal*csal)  +  (hLplus  (2 , 2 ) *snalS) ) ; 

PlusMet2=  . 5* ( (hLplus (1, 1) *csa2S)  +  (hLplus (1 , 2 ) *csa2*sna2 )  +  ... 

(hLplus (2, 1) *sna2*csa2)  +  (hLplus (2 , 2 ) * sna2S ) ) ; 

PlusMet3=  . 5* ( (hLplus (1, 1) *csa3S)  +  (hLplus ( 1 , 2 ) *csa3*sna3 )  +  ... 

(hLplus (2, 1) *sna3*csa3)  +  (hLplus (2 , 2 ) * sna3S ) ) ; 

ScalarMetl=  . 5* ( (hLscalar ( 1 , 1 ) *csalS )  +  (hLscalar ( 1 , 2 ) *csal* snal )  +  ... 

(hLscalar (2, 1) *snal*csal)  +  (hLscalar (2 , 2 ) *snalS) ) ; 

ScalarMet2=  . 5* ( (hLscalar ( 1 , 1 ) *csa2S )  +  (hLscalar ( 1 , 2 ) *csa2 *sna2 )  +  ... 

(hLscalar (2, 1) *sna2*csa2)  +  (hLscalar (2 , 2 )* sna2S )) ; 

ScalarMet3=  . 5* ( (hLscalar ( 1 , 1 ) *csa3S )  +  (hLscalar ( 1 , 2 ) *csa3*sna3 )  +  ... 
(hLscalar (2, 1) *sna3*csa3)  +  (hLscalar (2 , 2 )* sna3S )) ; 

%  Resetting  variables  used  in  calculating  the  response 
%  function  to  zero  for  the  next  iteration. 

DXl sum=0 ; 

DX2  sum=0 ; 

DX3sum=0; 

DX12sum=0 ; 

DX123Msum=0; 

DX123Ssum=0; 

DXScalarl sum=0 ; 

DXScalar2  sum=0 ; 


DXScalar3sum=0; 

DXScalarl2sum=0 ; 

DXScalarM12  3sum=0 ; 

DXScalarS12  3sum=0 ; 

DXC  omb  o 1 s  um= 0 ; 

DXC  omb  o  2  s  um= 0 ; 

DXC  omb  o3sum=0; 

DXCombol2sum=0; 

DXComboM123sum=0; 

DXComboS123sum=0; 

%  Averaging  over  epsilon  accounts  for  rotation  about 
%  the  third  Euler  angle  required  in  transforming  from 
%  the  source  frame  to  the  heliocentric  frame, 
for  epsilon  =  epsilonl 

%  Averaging  over  epsilon  accounts  for  the  rotation  about  the  third 
%  Euler  angle. 

metepsone  =  cos (2*epsilon) *CrossMetl  -  sin ( 2*epsilon) *PlusMetl ; 
metepstwo  =  cos (2*epsilon) *CrossMet2  -  sin ( 2*epsilon) *PlusMet2 ; 
metepsthree  =  cos (2*epsilon) *CrossMet3  -  sin (2*epsilon) *PlusMet3; 

ComboCrossScalarl  =  metepsone  +  Mixf actor *ScalarMetl ; 

ComboCrossScalar2  =  metepstwo  +  Mixf actor *ScalarMet2 ; 

ComboCrossScalar3  =  metepsthree  +  Mixf actor *ScalarMet3 ; 

%  Resetting  values  used  in  delta  loop. 

DX1=0; 

DX2=0; 

DX3=0; 

DX12=0; 

DX123M=0; 

DX123S=0; 

DXScalarl=0 ; 

DXScalar2=0; 

DXScalar3=0 ; 

DXScalarl2=0; 

DXScalarM123=0; 

DXScalarS123=0; 

DXCombol=0; 

DXCombo2=0; 

DXCombo3=0 ; 

DXCombol2=0; 

DXComboM123=0; 

DXComboS123=0; 

%  Delta  is  the  phase  constant  for  the  incident  wave,  (kx-wt+delta) 

%  This  portion  of  the  program  finds  the  amplitude  of  the  response. 

%  The  response  will  vary  depending  upon  the  initial  phase  of  the 
%  gravitational  wave.  As  the  program  cycles  through  possible  values 
%  for  the  initial  phase,  it  stores  the  maximum  response, 
for  delta  =  pi/32 : pi/ 1 6 : 31*pi/32 ; 

%  fl,  f2,  and  f3  include  the  time  dependent  parts  for  the  response 
%  function. 

fl=  ( (sin (delta-muone*u) -sin (delta-u) ) / (u* (1-muone) ) + (sin (u+delta) 
-sin (delta-muone*u) ) / (u* ( 1+muone) ) ) ; 
f2=  ( (sin (delta-mutwo*u) -sin (delta-u) ) / (u* (1-mutwo) ) + (sin (u+delta) 
-sin (delta-mu two* u) ) / (u* ( 1+mutwo) ) ) ; 
f3=  (  (sin (u* (1+muone) -delta) -sin (mutwo*u-delta) ) / (u* (1-muthree) )  . . 


+ (sin (u* (1-muone) +delta) +sin (mutwo*u-delta) ) / (u* (1+muthree) ) ) ; 


%  deltaX  values  correspond  to  the  change  in  length  of 
%  each  arm.  This  combines  the  time  independent  components 
%  with  the  time  dependent  component  to  achieve  the  complete 
%  response  of  LISA  to  the  gravitational  wave. 

%  XI  corresponds  to  LI,  X2  to  L2,  and  so  on. 

%  Similar  calculations  were  computed  for  the 
%  scalar  response  and  the  response  to  the 
%  combination  of  scalar  and  GR  polarizations. 
deltaXl=  metepsone* f 1 ; 
deltaX2=  metepstwo*f2 ; 
deltaX3=  metepsthree*f3; 

%  In  order  to  calculate  the  response  we  must  find  the  delta  that 
%  yields  the  maximum  result  (indicating  amplitude  of  response) . 
if  abs(deltaXl)  >  DXl 
DXl  =  abs (deltaXl ) ; 

end 

if  abs(deltaX2)  >  DX2 
DX2  =  abs (deltaX2 ) ; 

end 

if  abs(deltaX3)  >  DX3 
DX3  =  abs (deltaX3 )  ; 

end 

if  abs (deltaXl-deltaX2 )  >  DX12 
DX12  =  abs (deltaXl-deltaX2 ) ; 

end 

if  abs (deltaXl+deltaX2-2*deltaX3)  >  DX123M 
DX123M  =  abs (deltaXl+deltaX2-2*deltaX3) ; 

end 

if  abs (deltaXl+deltaX2+deltaX3)  >  DX123S 
DX123S  =  abs (deltaXl+deltaX2+deltaX3) ; 

end 

deltaXScalarl  =  ScalarMetl*fl ; 
deltaXScalar2  =  ScalarMet2*f2 ; 
deltaXScalar3  =  ScalarMet3*f3; 

if  abs (deltaXScalarl )  >  DXScalarl 
DXScalarl  =  abs (deltaXScalarl ) ; 

end 

if  abs (deltaXScalar2 )  >  DXScalar2 
DXScalar2  =  abs (deltaXScalar2 ) ; 

end 

if  abs (deltaXScalarS)  >  DXScalar3 
DXScalar3  =  abs (deltaXScalarS ) ; 

end 

if  abs (deltaXScalarl-deltaXScalar2 ) >  DXScalarl2 

DXScalarl2  =  abs (deltaXScalarl-deltaXScalar2 ) ; 

end 

if  abs  (deltaXScalarl+deltaXScalar2-2 *deltaXScalar3 )  >  DXScalarM123 
DXScalarM123  =  abs  (deltaXScalarl+deltaXScalar2-2 *deltaXScalar3 ) 

end 

if  abs  (deltaXScalarl+deltaXScalar2+deltaXScalar3 )  >  DXScalarS123 
DXScalarS123  =  abs  (deltaXScalarl+deltaXScalar2+deltaXScalar3 ) ; 

end 

deltaXCombol  =  ComboCrossScalarl*fl ; 
deltaXCombo2  =  ComboCrossScalar2*f2 ; 
deltaXCombo3  =  ComboCrossScalar3*f3; 


if  abs (deltaXCombol )  >  DXCombol 


DXCombol  =  abs (deltaXCombol ) ; 

end 

if  abs (deltaXCombo2 )  >  DXCombo2 
DXCombo2  =  abs (deltaXCombo2 )  ; 

end 

if  abs (deltaXComboS )  >  DXComboS 
DXComboS  =  abs (deltaXComboS ) ; 

end 

if  abs (deltaXCombol-deltaXCombo2 ) >  DXCombol2 

DXCombol2  =  abs (deltaXCoinbol-deltaXCombo2 ) ; 

end 

if  abs  (deltaXCombol+deltaXCombo2-2*deltaXCombo3)  >  DXComboM123 
DXComboM123  =  abs  (deltaXCombol+deltaXCombo2-2*deltaXCombo3 ) 

end 

if  abs  (deltaXCombol+deltaXCombo2+deltaXCombo3)  >  DXComboS123 
DXComboS123  =  abs  (deltaXCombol+deltaXCombo2+deltaXCombo3 ) ; 

end 

end  %  End  delta  loop. 


%  Add  the  responses  for  each  arm  and  each  combination  of  arms. 
DXlsum  =  DXlsum  +  abs(DXl); 

DX2sum  =  DX2sum  +  abs(DX2); 

DX3sum  =  DX3sum  +  abs(DX3); 

DX12sum  =  DX12sum  +  abs(DX12); 

DX123Msum  =  DX123Msum  +  abs (DX123M) ; 

DX123Ssum  =  DX123Ssum  +  abs (DX123S) ; 

DXScalarlsum=DXScalarlsum  +  abs (DXScalarl )  ; 
DXScalar2sum=DXScalar2sum  +  abs (DXScalar2 ) ; 
DXScalar3sum=DXScalar3sum  +  abs (DXScalar3 ) ; 
DXScalarl2sum=DXScalarl2sum  +  DXScalarl2; 
DXScalarM123sum=DXScalarM123sum  +  DXScalarM123; 
DXScalarS123sum=DXScalarS123sum  +  DXScalarS123; 


DXCombolsum=DXCombolsum  +  abs (DXCombol ) ; 
DXCombo2sum=DXCombo2sum  +  abs (DXCombo2 ) ; 
DXCombo3sum=DXCombo3sum  +  abs  (DXCombo3) ; 
DXCombol2sum=DXCombol2sum  +  DXCombol2; 
DXComboM123sum=DXComboM123sum  +  DXComboM123 ; 
DXComboS123sum=DXComboS123sum  +  DXComboS123 ; 


end  %  End  epsilon  loop 

%  Stores  the  sum  of  the  responses  for  each  angle  sampled. 
A(N,n)  =  DXlsum; 

B  (N, n)  =  DX2sum; 

C  (N, n )  =  DX3sum; 

D(N,n)  =  DX12sum; 

E(N,n)  =  DX123Msum; 

F(N,n)  =  DX123Ssum; 


ScalarA (N, n) 
ScalarB  (N, n) 
ScalarC (N, n) 
ScalarD  (N, n) 
ScalarE  (N, n) 
ScalarF  (N, n) 


DXScalarlsum; 

DXScalar2sum; 

DXScalarlsum; 

DXScalarl2sum; 

DXScalarM123sum; 

DXScalarS123sum; 


ComboA(N,n)  =  DXCombolsum; 
ComboB(N,n)  =  DXCombo2sum; 
ComboC(N,n)  =  DXCombolsum; 
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ComboD(N,n)  =  DXCombol2 sum; 

ComboE(N,n)  =  DXComboM123sum; 

ComboF(N,n)  =  DXComboS123sum; 

end  %end  t  loop 

end  %end  u  loop 

%  Divide  by  the  number  of  samples  summed  to  achieve  an  averaged  GR 
%  polarization  (by  averaging  over  epsilon) . 

A  =  (1/ceps ) . *A; 

B  =  (1/ceps) . *B; 

C  =  (1/ceps) . *C; 

D  =  ( 1/ceps )  . *D; 

E  =  (1/ceps) . *E; 

F  =  (1/ceps ) . *F; 

ScalarA  =  ( 1/ceps ). *ScalarA; 

ScalarB  =  ( 1/ceps ). *ScalarB; 

ScalarC  =  ( 1/ceps ). *ScalarC ; 

ScalarD  =  ( 1/ceps ). *ScalarD; 

ScalarE  =  ( 1/ceps ). *ScalarE ; 

ScalarF  =  ( 1/ceps ). *ScalarF; 

ComboA  =  (1/ceps) . *ComboA; 

ComboB  =  (1/ceps) . *ComboB; 

ComboC  =  (1/ceps) . *ComboC; 

ComboD  =  (1/ceps) . *ComboD; 

ComboE  =  (1/ceps) . *ComboE; 

ComboF  =  (1/ceps) . *ComboF; 

%  Preparing  to  generate  a  loglog  plot 

bA  =  logic  (A) ; 

bB  =  logic  (B)  ; 

bC  =  logic  (C); 

bD  =  logic  (D) ; 

bE  =  logic  (E); 

bF  =  logic  (F) ; 

bSA=loglC  (ScalarA) ; 
bSB=loglC  (ScalarB) ; 
bSC=loglC  (ScalarC) ; 
bSD=loglC  (ScalarD) ; 
bSE=loglC  (ScalarE) ; 
bSF=loglC  (ScalarF) ; 

bCA  =  logic  (ComboA) ; 
bCB  =  logic  (ComboB) ; 
bCC  =  logic  (ComboC) ; 
bCD  =  logic  (ComboD) ; 
bCE  =  logic  (ComboE) ; 
bCF  =  logic  (ComboF) ; 

%  Generate  a  plot  to  signify  the  program  has  finished  running  and  allow  the 
%  user  to  verify  an  initial  set  of  the  results, 
figure 
surf  (bA) 

%  End  program  run  time  and  display  in  the  command  prompt, 
toe 

%  Saves  variables  for  reference  eliminating  the  need  to  run  the  program  more 
%  then  once  for  a  set  of  values, 
save  (  ' ThreeDplots31MAR ' ) ; 
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